ROME1/1436/06 



Modelling non-perturbative corrections 
to bottom-quark fragmentation 



Ugo Aglietti 12 , Gennaro Corcella 1 
and Giancarlo Ferrera 1 2 

1 Dipartimento di Fisica, Universita di Roma 'La Sapienza' 
P.le A. Mow 2, 1-00185 Roma, Italy 

2 IN FN, Sezione di Roma 
P.le A. Mow 2, 1-00185 Roma, Italy 



Abstract 

We describe i?-hadron production in e + e~ annihilation at the Z° pole by means of a 
model including non-perturbative corrections to 6-quark fragmentation as originating, via 
multiple soft emissions, from an effective QCD coupling constant, which does not exhibit 
the Landau pole any longer and includes absorptive effects due to parton branching. We 
work in the framework of perturbative fragmentation functions at NLO, with NLL DGLAP 
evolution and NNLL large- x resummation in both coefficient function and initial condition 
of the perturbative fragmentation function. We include hadronization corrections via the 
effective coupling constant in the NNLO approximation and do not add any further non- 
perturbative fragmentation function. As part of our model, we perform the Mellin transforms 
of our resummed expressions exactly. We present results on the energy distribution of b- 
flavoured hadrons, which we compare with LEP and SLD data, in both x- and A-spaces. 
We find that, within the theoretical uncertainties on our calculation, our model is able to 
reasonably reproduce the data at x < 0.92 and the first five moments of the B cross section. 



1 Introduction 



We study 6-quark fragmentation and 5-hadron production in e + e~ annihilation using a model 
which implements power corrections via an analytic effective coupling constant. Such a model will 
allow us to make hadron-level predictions, directly comparable with the accurate experimental data 
from the SLD [1] and LEP [2-4] collaborations. 

Let us describe in physical terms the process we are interested in: 

e+e" -> Z° -> B + X- b , (1.1) 

where X b is a generic hadronic final state with a b quark. The corresponding parton-level process, 
assuming that an arbitrary number of gluons is radiated, reads: 

Z° -> bb^bb 9l ...g n . (1.2) 

In the Z° rest frame, the b and the b are initially emitted back-to-back with a large virtuality, of 
the order of the hard scale Q = mz- Then, they reduce their energy by emitting gluons, mainly of 
low energy or at small angle (soft or collinear radiation). As long as the virtuality of the b quarks is 
large enough, we can neglect the b and b masses. Considering, for simplicity, the case of one single 
emission, 

e+ e - ^ Z \q) - b{p b ) + b(p- b ) + g(p g ), (1.3) 
and defining the energy fractions 



Xi = — = , 1.4 

m z mz 

where i = b,b,g and x b + x b + x g = 2, a straightforward kinematical computation yields: 

1-^ = ^(1 - 008^), (1.5) 

with dig being the angle between the gluon and the b. From Eq. (11. 5p . we see that collinear emission 
from the b (9i g = 0) implies Xb=l- Likewise, by crossing b —>b, one can show that collinear radiation 
off the b, i.e. d bg = 0, corresponds to x b = 1. Soft-gluon radiation {x g = 0) yields x b = xi = 1. We 
can identify two main different mechanisms of energy loss of the b: a direct loss, related to soft or 
collinear emissions off the b itself, and an indirect loss, when the gluons are radiated off the b (hard 
and large-angle emissions are not enhanced). The asymmetry between the b and the b jets is not 
dynamical, but is just related to the fact that, e.g., we decided to measure the energy of the b, and 
not that of the b. 

When the virtuality of the b has reduced to a value of the order of its mass, mass effects become 
substantial, and when it becomes comparable with the hadronic scale, 0(1 GeV), the b approaches 
the non-perturbative phase of its evolution. The hadronization, for example, into a B meson can 
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be described as due to the radiation of a light qq pair, and the subsequent combination of the b 
and the q to form a B = (bq) state: 

b -> B + q. (1.6) 

In the process (11 .6ft . the colour of the heavy 6 is coherently transferred to the light g; we shall assume 
that no more radiation is emitted after the colourless hadron has been formed. The hadronization 
transition (11.61) is sensitive to non-perturbative power corrections A/m&, with A being the QCD 
scale, e.g., in the MS scheme. Such power corrections - for the time being - cannot be calculated 
from first principles, but are usually described by means of phenomenological hadronization models, 
such as the Kartvelishvili [5] or the Peterson [6] models, containing parameters which are to be fitted 
to experimental data. In Ref. [7], furthermore, non-perturbative effects in bottom fragmentation 
at large x were described as a shape function of m&(l — x) and implemented using Dressed Gluon 
Exponentiation [8]. 

In this paper, we shall follow a different approach to model power corrections to 6-quark frag- 
mentation. We shall assume that the transition (11.61) can be described in terms of an effective 
analytic coupling constant as, which incorporates non-perturbative power corrections and, unlike 
the models above mentioned, does not present any free parameter to be fitted to data. Within our 
model, whenever we use as instead of the standard as, the energy of the b quark will be identified 
with the one of the observed B meson, i.e. E\, —>■ Eg. The non-perturbative model which we 
shall use hereafter was already used in [9] in the context of 5-meson decays and interesting results 
were found. The comparison with the data on the photon spectrum in radiative decays and on the 
hadron- mass distribution in semileptonic decays showed good agreement, while discrepancies were 
found for the purpose of the electron energy distribution. 

An analogous model can be introduced for the fragmentation of a b into a baryon, such as, for 
example, a hyperion A&. The b quark emits two light uu and dd pairs and combines with the u and 
the d to form the colourless At = (bud): 

b -y A b + u + d, (1.7) 

where the ud system is emitted with the same colour state as the initial b. The process (11.71) will 
be also described by means of an effective coupling constant: however, there is no reason a priori 
to assume that the effective coupling constant for (jl.7p is equal to that for (II. 6p . 

Turning back to the parton-level process (ll.2p . the differential cross section has a perturbative 
expansion containing large logarithms of the form: 

2 

777. 

a n s ln fc — § (n = 1, 2, 3 • • • oo, k = 1, 2, • • • , n), (1.8) 

m b 

which are related to collinear emissions of partons with transverse momenta ranging between scales 
fixed by the particle masses: 

. ml f m l dk\ .„ _ 

a s m— f = a s I 1.9 
mi Jm\ k L 
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At higher orders, the leading contributions, k = n in Eq. (11.81) . are associated with phase-space 
regions ordered in transverse momentum: 



> fcl2 > 



leading to integrals of the form 

a 



z dk 2 ±1 f k li dk 2 ±2 



2 h 2 
m b K ±l "' 



2 P 
b K ±2 



> k±nj 



i dk^ n 



lag ln"^f. 



1.10) 



1.11) 



Hence, in the first-order cross section, the collinear logarithm exponentiates. 
The 6-quark spectrum is also affected by large logarithms 

(n = l,2---,oo, fc = 0,l,2, 



ln fe (l -x) 



1 — x 



2n-l), (1.12) 



also called threshold logarithms, which are enhanced for x — > 1, i.e. for soft or collinear emission. 
The plus distributions in Eq. fll . 12j) are defined as: 



"ln fc (l -x) 


= lim 


1 — X 





9(1 -x-e 



\n k (l-x) 



x 



- 6(1 -x -e) 



i-« ln fc (l -x') 



fix' 
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The contributions in (I1.12p originate from the following double-logarithmic integrals: 

duo dt 

o(l — x — a;£J 



ln(l — x) 




17 


1 — X 


+ 


Jo JO 



(1.13) 



1.14) 



where cj = 2E g /mz, is the normalized energy of the radiated gluon and t = 2(1 — cos 6) ~ with 
9 being the emission angle. The 5-function contains a typical kinematical constraint and the 'plus' 
regularization comes after including the virtual diagrams. To factorize the kinematical constraint 
for multiple gluon emissions, a transformation to iV-space is usually made: 







Cdxx N - 1 


/7 


h 


JO Jo 



5(1 — x — ut) 



1 r 1 duo dt 



o Jo to t 



(1 - cot) 



N-l 



~ - - In 2 N. 



(1.15) 

At higher orders, the large-x logarithms are associated with multiple soft emissions ordered in 
angle: 

0x > 2 > ••• > e n . (i.i6) 

Colour coherence [10], in fact, dictates that multiple soft radiation interferes destructively outside 
the angular-ordered region (11.161) of the phase space. That produces integrals in moment-space of 
the form 



a c 



n f 1 f 1 dtoi dt\ f 1 f 1 du} 2 dt 



JO LOi t\ JO JO L0 2 t 2 
N-l 



1 r 1 doo„ dt r 



X 



1 - uih) 

( r 1 r 1 duo dt 
Wo Jo CO t 



1 



JO LO n t n 

N-l 



e(ti > t 2 > ■ ■■ > t n ) 



00 t) 



1 - L0 2 t 2 
N-l 



- 1 



C^n tn 



N-l 



1.17) 



with uji = lEijmz and and U = (1 — cos6>j)/2 ~ Of. Eq. (11.171) implies that also the threshold 
logarithms have an exponential structure. Schematically, the inclusion of higher orders amounts to 
the replacement: 

1 + a s In 2 AT -> e asln2N . (1.19) 

A consistent method to accomplish the resummation of mass (ll.§p and threshold (11.121) logarithms 
is the formalism of perturbative fragmentation functions [11]. The basic idea is that a heavy quark 
is first produced at large transverse momentum, k± ^> m, and can be treated in a massless fashion, 
afterwards it slows down and fragments into a massive parton. This leads to writing the energy 
distribution of the b quark, up to power corrections, as the following convolution: 

~ m 2 z , m 2 b ) = C(x; m%, /4) ® D(x; /j? f , m 2 b ), (1.20) 

where: 

1. C(x; m|, is a coefficient function, obtained from a massless computation in a given fac- 
torization scheme, describing the emission off a light parton. The coefficient function contains 
large-x logarithms, as in Eq. (11.121) . which are process-dependent; 

2. D(x; fi F ,m 2 ) is the perturbative fragmentation function, associated with the transition of a 
massless parton into the heavy b. 

The perturbative fragmentation function follows the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi 
(DGLAP) evolution equations [12,13], which can be solved once an initial condition is provided. 
Solving the DGLAP equations one resums the large logarithms (jl.8p which appear in the massive 
6-quark spectrum. Function D(x; fi F ,rrib) can therefore be factorized as: 

D{x- /4, m 2 b ) = E{x- /4, fi 2 QF ) ® D ini (x; n 2 QF , m 2 b ), (1.21) 

where: 

1. E{x; fi 2 F , Hq F ) is an evolution operator from the scale fj, F , typically of the order of the hard 
scale mz, down to [Iq F , of the order of the bottom-quark mass m^. In E(x; fi F , fJ^ F ), the mass 
logarithms (11. 8ft are resummed; 

1 The integral on the l.h.s. of Eq. (|1.17j) is completely factorized into single-particle integrals except for the 
angular-ordering 0-function. This constraint can be eliminated by symmetrizing over the t^s: 

e(ti > h > ■■■ > i n ) -». — e (^l > *<a > ••• > tin) = -■ (1-18) 

perm 
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2. D im (x; fiQ F , mf ) is the initial condition of the perturbative fragmentation function at the scale 
Hqf — ma, first calculated in [11], and lately proved to be process-independent [14]. It also 
contains threshold logarithms fll.l2p . which are also process- independent, and were resummed 
in [14] in the NLL approximation. 

The approach of perturbative fragmentation functions, with the inclusion of NLL large-x re- 
summation, has been applied to investigate 6-quark production in e + e~ annihilation [14], top 
(t — > bW) [15,16] and Higgs decays (H — > bb) [17]. In this paper we shall go beyond the NLL 
approximation and, as far as the perturbative calculation is concerned, we shall also include next- 
to-next-to leading logarithmic (NNLL) threshold contributions to both e + e~ coefficient function 
and initial condition of the perturbative fragmentation function. Also, the Mellin transforms of our 
resummed expressions will be performed exactly, in order to resum constants and power-suppressed 
terms as well. 

Most previous analyses convoluted the parton-level spectrum with a non-perturbative frag- 
mentation function, which contains few parameters which are to be fitted to experimental data. 
Afterwards, such models are used to predict the B spectrum in other processes, as long as one 
consistently describes perturbative b production (see, e.g., Ref. [18]). An alternative method [19] 
consists in using data in Mellin moment space, such as the DELPHI ones [4], and fit directly the 
moments of the non-perturbative fragmentation function, without assuming any functional form 
for the hadronization model in x-space. 

In this paper we reconsider the inclusion of non-perturbative corrections to bottom-quark frag- 
mentation. As discussed before, instead of fitting the parameters of a hadronization model, or the 
moments of the non-perturbative fragmentation function, we model non-perturbative effects by the 
use of an analytic effective coupling constant [20-22], which does not contain any free parameter, 
extending the analysis carried out in [23] and in [9] in the framework of heavy-flavour decays. Our 
modified coupling constant does not exhibit the Landau pole any longer, and includes all-order ab- 
sorptive effects due to gluon branching. We shall then be able to compare our predictions directly 
with data, without using any extra hadronization model. 

The plan of our paper is the following. In Section [2] we review the main features of the MS 
e + e~ coefficient function. In Section [3] we discuss the massive computation, the perturbative 
fragmentation approach and the resummation of the large mass logarithms. In Sections H] and 
we discuss NNLL large-x resummation in the coefficient function and initial condition of the 
perturbative fragmentation function, respectively, pointing out the differences of our analysis with 
respect to previous ones. In Section [6] we construct the analytic coupling constant without the 
Landau pole. In Section [7] we present our model, based on an effective coupling constant as the 
only source of non-perturbative corrections to 6-quark fragmentation. In Section [8] we present our 
results on the 75-hadron spectrum in e + e~ annihilation in x-space, and investigate how they fare 
against SLD, ALEPH and OPAL data. In Section [9] we perform a similar analysis in iV-space and 
compare with the DELPHI moments. In Section [10] we summarize the main results of our work and 
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discuss the lines of development of our study. There is also an appendix describing an algorithm 
to compute numerically the Mellin transforms with the Fast Fourier Transform. 



2 Massless- quark production and NLO coefficient function 



In this section we consider massless-quark production and discuss the main features of the NLO 
e + e~ coefficient function, which will be used later on in the framework of perturbative fragmentation 
functions. The NLO computation exhibits many properties of the higher orders, so that the general 
structure of the perturbative corrections can be understood by looking in detail into this case. 

We study the production of a light-quark pair at the Z° pole, in the NLO approximation: 

e+e- - Z°(q) - q(p q ) + q{p q ) + (g( Pg )), (2.1) 

where (g(p g )) denotes a real (virtual) gluon, and define the light-quark energy fraction as: 

2p 9 • q 



x = 



Tflr 



(2.2) 



The differential cross section, in dimensional regularization, can be read from the formulas in [24]: 



1 da i _ 5 n _ x \ + ?s(jj>r) C f 



a dx 



7T 





"1 + x 2 ~ 




1 — X 



m 



ln^f - -\+A(x)\+0(a 2 ), (2.3) 

f^F e / ) 



where 



11 , , A s 

- = 7p + ln(47rj 

e e 



(2.4) 



e = (4 — D)/2, D is the number of space-time dimensions, fin and /xp are the renormalization 
and factorization scales, je = 0.577216 ... is the Euler constant, as is the dimensionless strong 
coupling constant, Cf = 4/3 and A(x) is the following function, independent of e and /ip: 



A{x) 



ln(l — x) 
1 — x 



71 



J + 



4 [1 



x\ 



a 



T -3 *(1 



x) — — ln(l 



, 1 + x 2 , 5 - 3x 
x) H In x + 



1 — x 



00 



l 3 as(n 2 R )C F 

4 7T 



0(a 



4 
(2.5) 

(2.6) 



is the NLO total cross section, with <Jo being the Born one. Eq. (I2.3P presents a pole, ~ 1/e, which 
is remnant of the collinear singularity, and disappears in the total cross section. Subtracting the 
collinear pole, one obtains what, in the perturbative fragmentation formalism, is called the MS 
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coefficient function [11]: 



1 da„ 



-i MS 



5(1 -x) + 



a dx 

The coefficient function presents: 



<xs(t*k)C F Jl 


1 + x 2 ' 


71 [2 


1 — X 



rrir 



hx'-^f + A(x)\ + 0(a 2 s ). (2.7) 
\i F J 



1. a term of collinear origin, 



/4 



2tt 



1 + a; 2 
1 — x 



(2.8) 



where (x) is the leading-order Altarelli-Parisi splitting function, containing contributions 
enhanced in the threshold region x — > 1; 



2. two terms which are enhanced at large x and are independent of fi F : 



7T 



ln(l — x) 

1 — x 



and 



3 C F 1 

4 7T [1 - x] 



(2.9) 



3. a term proportional to 8(1 — x), i.e. a spike in the elastic point x — 1: 

— y-3U(l-x); 



(2.10) 



4. terms dependent on x and divergent at most logarithmically for x — > 1: 

a s C F 



7T 



1 + x 1 + x 2 5 — 3x 

— ^— ln(l - x) + y— - lnx + — ^— 



(2.11) 



We observe that the latter contribution also contains a term, ~ lnx, which is enhanced at 
small x. 



We wish to rearrange the coefficient function and put it into a form which will become more 
useful in the following sections, when dealing with large- x resummation. To this goal, we rearrange 
the Altarelli-Parisi splitting function by using the identity: 



1 + x 2 
1 — x 



[1 



x\ 



[l + x) + -6(1 -x). 



(2.12) 



2 Unlike the usual convention, where the coefficient function is expressed in terms of 1 / <Jo(da / dx) , we have 
divided by the NLO cross section a, so that the total integral is 1. In the following, we shall compare with data, 
also normalized to 1. 
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Moreover, we factorize the constants which multiply the term 
coefficient function in the form: 



8(1 — x), and finally write the MS 



' 1 da q 
a dx 



MS 



1 + 



a s (fJL 2 R )C F ( tc 2 



71 



4 n 2 Fi 



8(1 -x) + 



71 



4(1 



x 



a s (fx 2 R )C F jl + x 2 



71 



1 — X 



Inx 



1 + x 



ln(l 



x) 



In 



i4 



+ 



1 — X 
5 — 3x 



In 



m|(l 



x 



+ 0(a 2 s ), 



which is equal to (12.131) at 0(ots). Eq. ( 12. 13ft exhibits the logarithm ln[m|(l 
originates from the integration over the variable k 2 = (p q +p g ) 2 (l 



x) 



a, ln m z (l-x) 



7T 



Hp 



k 2 



71 



(2.13) 
x)//i F ], which 



(2.14) 



where the constant term Cpas/it has been brought inside the integral for reasons which will become 
clear later. For soft and collinear radiation, k 2 ~ E 2 8 2 g = k\, the gluon transverse momentum 
with respect to q. In principle, the lower value of k 2 would be zero, as massless quarks can emit 
soft gluons at arbitrarily small angles. In dimensional regularization, using the MS factorization 
scheme, the minimum k 2 is set by the factorization scale: k^ in = fi 2 F . The upper limit in the integral 



(I2.14p can be obtained observing that k 2 = (p q + p g ) 



x) 



(q-P<i) (1 -x)< m%(l - x) 



3 Heavy-quark production 

and NLO perturbative fragmentation function 



Let us now consider the production of massive bottom quarks in the NLO approximation: 

> Z°(q) - b(p b ) + l( Pl ) + (g(p g )). (3.1) 



e e 



The differential cross section for the production of a b quark of energy fraction x reads [11]: 
1 dab 



a dx 





"1 +x 2 ~ 


In— | - 

mi 


ln(l — x) 


7 1 


7T \2 


I — X 


1 — X 


+ 4(l-x) + 



+ 



8(1 -x) + 

( k 2 \ „, . 1 + x , , 
__ 2 8(l-x) + —\a(l-x) 



(3.2) 



1 



x 



Inx + 



7 — x 



+ o 



x 



®S(H 2 R, 



where p > 1. The massive spectrum, unlike the massless one, is free from collinear singularities 
because the quark mass acts a regulator. However, Eq. (13.21) presents a large mass logarithm, 
~ as ln(m%/m 2 ), which needs to be resummed to improve the perturbative prediction. 



S 



The resummation of ln(m 2 z /ml) can be achieved by the use of the approach of heavy-quark 
perturbative fragmentation functions [11], which factorizes the spectrum of a massive quark as the 
following convolution: 



1 da, 
a dx 



b i 2 2^ 

a;; m z ,m b ) 



E 



1 dz 

x Z 



1 d&i 
a dz 



(z,m|,/4) 



MS 



D 



MS 



x 



; fjL 2 F ,m 2 b ) + O 



fm h \ p ' 
\m z 



(3.3) 



In Eq. (13.31) . 1/er (d&i/dx) is the coefficient function, corresponding to the production of a massless 
parton i, Di(x,m b , fi F ) is the heavy-quark perturbative fragmentation function, associated with 
the transition of the massless parton i into a heavy b, and fi F is the factorization scale. In the 
following, we shall neglect b production via gluon splitting g — > bb, which is negligible at the Z° 
peak and suppressed at large x. Hence, in Eq. ( 13 . 3[) . i = b and 1/a ida b jdz) is the quark coefficient 
function, presented in Eq. (12. 7p in the MS factorization scheme. The perturbative fragmentation 
function _D^ S expresses the fragmentation of a massless b into a massive b. 

Requiring the massive cross section to be independent of fi F , one obtains that the perturbative 
fragmentation function follows the DGLAP evolution equations [12, 13], which can be solved once 
an initial condition is given. The initial condition of the perturbative fragmentation function 
D b m (x b , Hq F , m 2 ) was given in [11], and can be obtained inserting in Eq. (13.31) the massive spectrum 
(EQ) and the MS coefficient function ffOh It reads: 



D™(x; a s (fi 2 0R ),fi 2 0R ,fi 2 0F ,m 2 b ) = 5(l-x)4 



71 



1+x 2 
1 — X 



5*4- Mi-*) 



+ 0{a 2 s ). 
(3.4) 

In [14], the process- independence of the initial condition (13.41) was established on general grounds. 
The solution of the DGLAP equations is typically obtained in Mellin moment space. At NLO, and 
for an evolution from /x f to fi F , it is given by: 



D KN (fx 2 F ,m 2 b ) = E N a s {iil F ),as{ii F ) 



u b,N 



<Xs((4r),(J%r,(J>If,™. 



(3.5) 



where D b m N is the Mellin transform of Eq. (13. 4p and En is the DGLAP evolution operator [12, 13]: 



E 



N 



(!/■ h ®s(Hf) 



(0) 



exp 



N 



hi 



If) 



+ 



(*s((J-of) ~ a s{^ F ) 
4tt 2 /5 



N 



2tt(3 



1 p(0) 
~ N 



(3.6) 



In Eq. (13. 6p . Pffl and are the Mellin transforms of the LO and NLO splitting functions; (3$ 
and /?i are the first two coefficients of the QCD /5-function: 

33 - 2n f 153 - 1971/ 



0o 



01 



12tt ' , ' ± 24tt 2 ' 

which enter in the NLO expression of the strong coupling constant: 

1 f 1 In [ln(Q 2 /A 2 )] 



0oHQ 2 /A 2 ) 



2 o HQ 2 /A 2 ) 



(3.7) 



(3.8) 
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As discussed in [11], Eq. (13.51) resums leading (LL) ag(fi 2 F ) \ia n (fip / fi^p) and next-to-leading (NLL) 
ag(fip) In 11-1 (fip / nl F ) logarithms of the ratio of the two factorization scales (collinear resumma- 
tion). For fip ~ mz and //qf — m b, as we shall assume hereafter, one resums the large logarithm 
ln(m|/7??|), which appears in the NLO massive spectrum (13. 2p . with NLL accuracy. 

As we did for the coefficient function, we rearrange the initial condition of the perturbative frag- 
mentation function into a form which will be convenient when discussing soft-gluon resummation. 
We use the identity (I2.12p and the relation 



1 + x 2 



X 



ln(l 



x) 



= 2 


ln(l — x) 


+ 


1 — X 



J + 



1 + x) ln(l - x) - -5(1 



x\ 



(3.9) 



Furthermore, we factorize the coefficient of the term ~ 5(1 — x) and write D b m in the following 
form, which is equivalent to (13. 4p . up to terms of 0(a 2 s ): 



D b m [x;as(fil R )^l F: m 2 b } 



+ 



+ 



1 + a s(HoR)C F L 3 ^ nip ** 
it \ 4 m 2 b/ 



6(1 



®s(hIr)Cf 



In 



x 



ml(l — x) 2 



(1-x) 



®s(v r) c 



[1 + x) 



71 



2 2 ml V ; 



+ O(c4).(3.10) 



The logarithm \n[/iQp / (m b (l — x) 2 )] comes again from an integral over k 2 = (pb + p g ) 2 (l 

dk 2 



x) 



Up — m — ^ 

7r mi (1 — xr 



4f 



77 ^F • 

2 b {i-x) 2 k 71 



(3.11) 



As in (I2.14p . for soft and small-angle radiation k 2 ~ k\, the transverse momentum of the gluon 
relative to the b. Following [14] , the lower limit of the /c 2 -integration is easily found by considering 
the dead-cone effect [10]. In fact, unlike the coefficient function, where quarks are treated as 
massless, soft radiation off a massive b quark is suppressed at angles lower than [§: 



H ~ 

"mm — • 



(3.12) 



To the logarithmic accuracy we are interested in, we can neglect emission inside the dead cone and 
use Eq. (I3.12p to obtain the lower limit on the emitted-gluon transverse momentum: 



min J_min 



E 2 6 2 ■ ~ 

g "mm 



m 2 (l — x) 2 . 



(3.13) 



Comparing Eq. (13.111) with (12.141) . we observe that, in the MS scheme, the factorization scale 
squared fi 2 F is the lower limit on the gluon transverse momentum in the coefficient function; /Xq F is 



3 The bremsstrahlung spectrum off a heavy quark reads: da ~ as (dui/uj) [d6 2 /(0 2 + ^ ml )], where u> is the energy 
of the radiated soft gluon. 
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instead the upper limit on k 2 in the initial condition. The operator E(as(/J I Q F ), as(fJ, F )), given in 
Eq. (13. 6p describes the evolution between these two scales. 

Before closing this section, we point out that, although we shall use NLO coefficient function 
and initial condition, and evolve the perturbative fragmentation function with NLL accuracy, in 
principle we could go beyond such approximations. In fact, the NNLO e + e~ coefficient function was 
calculated in [25,26] and NNLO corrections to the initial condition of the perturbative fragmen- 
tation function in [27,28]. Furthermore, NNLO contributions to the non-singlet time-like splitting 
functions were computed in [29], while in the singlet sector they are still missing. In any case, we 
shall delay the inclusion of such NNLO corrections to future work. 



4 Large-ir resummation in the coefficient function 

In this section we perform threshold resummation in the e + e~ coefficient function to next-to-next- 
to-leading logarithmic accuracy, and combine the resummed result with the first-order one presented 
in Section [2j 



4.1 Resummed coefficient function in Mellin space 



According to Eq. (12.71) . the NLO e + e coefficient function contain terms of the form 



ln(l — x) 
1 — x 



and 



a s 



[l-x] 



(4.1) 



which become large in the limit x — > 1, corresponding to soft or collinear gluon radiation. All-order 
resummation in the coefficient function can be performed following the general lines of Refs. [30,31], 
and was implemented in [14] in the NLL approximation. 

Threshold resummation is typically performed in iV-space [14,31], where kinematical constraints 
factorize and the x — > 1 limit corresponds to iV — > oo. The Mellin transform of the NLO coefficient 
function, which we denote by Cjv, can be found in [11]. At large N, it exhibits single and double 
logarithms of the Mellin variable: 



C 



N 



1 + 



71 



1 1 2 *r ( 3 , mi 



/4 



In AT 



+ Q{^ F ,m 2 z ) + O 



N 



with the constant terms given by: 



Q(/4, m|) = ( ~ - 7B ) In ^ + ^vr 2 - 3 + ]p 2 E + ^ lE . 



(4.2) 



(4.3) 
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Besides the contributions (14.1j) . the NLO coefficient function (12.71) also contains a contribution 
~ ln(l — x). However, its Mellin transform behaves like ~ In N/N and is therefore 0(1/ N) at large 
N. The other terms in Eq. (12. 7p are suppressed at large N. 



The resummed coefficient function has the following generalized exponential structure: 



A 



(C) 
N 



/0\00 O C ( C~ 1 '\ / O \ o o o "1 



(4.4) 



where [14] 



G 



A' 



as(/4)>/4>*4>"4 



The exponent G^ 7 [asO"^); ^fi m %\ resums the large logarithms appearing in Eq. (14.21) : 



1 Z N-1 

dz 

o 1 — 



m%(l-z) dk 2 

m 2 f k' 



-A 



a s (k 2 ) +B 



a s 



— z 



n 1 n+1 



LL : a n s In 



N\ 



NLL : a n s In" JV; 



NNLL : al ln n_1 iV. 



(4.5) 



(4.6) 



As in [31], the integration variables are z — 1 — x g , x g being the gluon energy fraction, and 
k 2 = (Pb + Pg) 2 (l — z). In soft approximation, z ~ x; for small-angle radiation k 2 ~ k\, the gluon 
transverse momentum with respect to the h. 



The functions A(ots) and B(as) can be expanded as a series in as as: 



n=l 



n=l 



5(a5) = Effr^ ) 



(4.7) 
(4.8) 



In the NLL approximation, one needs to include the first two coefficients of A(as) and the first of 
B(as); to NNLL accuracy, A^ and are also needed. Their expressions read: 



A W 

A (2) 
A (3) 



C F 
1 



67 
18 



C(2) 



9 



£>(!) 
S (2) 



G F |g a 

C F nf 
3 



245 11 . , 67 >/ x 11 A/ x 
^ + 24 C ' 3 )-S C(2) + T C(4) 



G A n/ 



209 7 >/ , 5 A/ x 
1432 + I2 C < 3 > " I8 C(2) 



"55 


C(3)" 


njl 


96 ~ 


2 


108 J 



; V 432 6 



(4.9) 
(4.10) 



(4.11) 
(4.12) 

(4.13) 
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where Ca = 3, nj is the number of active flavours and C,{x) = Y^=\ ^/n x is the Riemann Zeta 
function. The first two coefficients of function A(as) have been known for long time [31]; more 
recent is the calculation of [32]. Function B(a>s) is associated with the radiation off the 
unobserved massless parton, e.g. the b if one detects the b: the coefficient B^ was given in [31], 
B^ was computed in [37]. 

We can already observe that the integral over the transverse momentum in Eq. (14.51) is a 
generalization of the NLO integral in Eq. (12.141) . where 



A 



TT 



TT 



TT 



a s (k 2 



Analogously, function B(as) generalizes the coefficient of the ~ 1/(1 — x) + term in (12. 13ft : 



3 -C F ^ = B^ 



B 



TT 



TT 



X 



(4.14) 



(4.15) 



A delicate point of our approach concerns the integral over z in Eq. (14. 5p . Calculations which use the 
standard coupling constant as{k 2 ), such as [14], typically perform the ^-integration approximating 
the term ~ (z N ~ l — 1) in such a way that only logarithmically-enhanced contributions ~ ln fc N are 
kept in the exponent. To NLL accuracy, such an approximation reads [31]: 



1-2- 



-7b 



N 



(4.16) 



where Q(x) is the Heaviside function. Beyond NLL, the prescription (14.161) can be generalized as 
discussed in [33]. In fact, the integrations in z and k 2 in (14.51) involve the infrared region z — > 1 for 
any value of N. As observed in [34] , performing such integrations exactly, when using truncated 
expressions for A(as) and B(as), will lead to a factorial divergence, corresponding to a power 
correction which, in our case, would be ~ A/m^. Ref. [35], however, remarkably pointed out 
that this contribution is spurious as it is actually related to the fact that one employed truncated 
expressions for A(a,s) and B(as)- It was shown in [35] that the higher-order coefficients of such 
functions also present factorial singularities leading to contributions ~ A/rriz, which cancel the 
analogous terms coming from the exact integration, so that one is just left with a power correction 
~ A 2 /m|. It was hence argued on general grounds in [33,36] that, when considering truncated 
expressions for functions A(as) and B(as), employing the step function (I4.16P or its extensions 
beyond NLL [33] is a better approximation than performing the z-integration exactly, since the 
transformation (I4.16P and its generalizations do not lead to any spurious factorial growth of the 
cross section. 

However, all such results are valid as long as one uses the standard coupling constant in the 
resummed formulas. As we shall detail in Sections [6] and 0, we will model non-perturbative correc- 
tions by means of an effective QCD coupling constant, based on an extension of the work in [20-22], 
which does not present the Landau pole any longer and includes power-suppressed contributions. 
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As a result, some power corrections will be unavoidably transferred to the resummed coefficient 
function when employing the effective coupling in expressions like Eq. (14. 5B . independently of how 
one performs the longitudinal-momentum integration. To our knowledge, there is currently no anal- 
ysis, such as the ones in [35,36], on factorial divergences and power corrections in the coefficient 
function (14. 5p . when using an effective analytic coupling constant along with truncated A(as) and 
B(as)- We cannot therefore draw any firm conclusion on whether, within our model, it is better to 
perform the ^-integration in an exact or approximated way. Given the phenomenological aim of the 
present paper, we prefer to postpone such an investigation. In the following, as already done in [9] 
when the effective coupling was used in the context of heavy-flavour decays, we shall present results 
obtained performing the Mellin transform (14.51) exactly. In other words, doing the z-integration in 
(14.51) exactly should be seen, for the time being, as part of the non-perturbative model which we 
shall propose. 

Due to its complexity, we cannot express the result of the exact longitudinal-momentum inte- 
gration (14 .5p in a closed analytic form, but we shall perform it numerically. The integration over k 2 
could in principle be made analytically, and the result expressed in terms of polylogarithms. For 
simplicity's sake, however, we perform also the /^-integration numerically. 

For < z < 1, the argument k 2 of as(k 2 ) in Eq. (14 .51) varies from zero to the hard scale squared 
m|; this implies that the number of active quark flavours does change in the fc 2 -integration. In the 
following, we shall include correctly the variation of rif at the quark-mass thresholds when doing 
the numerical integration. 

Let us now see with an explicit analytic computation the difference between an exact Mellin 
transform and an approximation using the step function. The expansion of the exponent of the 
resummed coefficient function reads, to 0(as)'- 





(4.19) 
(4.20) 



where 



d k+l lnT{x) 
dx k+1 



(4.21) 
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are the polygamma function and T(x) is the Euler Gamma function. Using the large- N expansions: 



Si (AO 
S 2 (N) 



In AT + lE + 0{ — 



71 



+ o 



(4.22) 
(4.23) 



we can write for large A^: 



A 



(CO 

N 



1 + 



«s(/4)cv 



as 



71 



-\n 2 N 
2 



+ Q'(ti 2 F ,m 2 z ) + 0(- 



N 



(4.24) 



where Q'{fi 2 F,m 2 z ) collects the constant terms 



^// 9 9n 71-2 7f 3 , m 2 7 

Q'^ 2 F ,m 2 z ) = — + -J + - lE - lEln-f. 



(4.25) 



By comparing Eq. (14.181) with Eq. (14.241) . we explicitly see that performing the Mellin transform 
exactly, rather than with the step-function approximation, amounts to including also constants 

(C) 

and contributions power-suppressed at large A" in the exponent G N . If one used the 6-function 
approximation, as in [14], the 0(as) expansion of the resummed coefficient function would con- 
tain only the logarithmically-enhanced contributions ~ as In A" and ~ as\n 2 N. The resummed 
coefficient function in the original x-space is finally obtained by an inverse Mellin transform: 



A (C) 



c+ioo 

c—ioo 2,7li 



as(^ R )^ 2 R,^F,m 2 z 



(4.26) 



where the (real) constant c is chosen so that all the singularities of A^P lie to the left of the 
integration contour. The inverse transform (I4.26[) is also made exactly, in a numerical way. 



4.2 Matching of resummed and NLO coefficient function 

We wish to implement the matching of the resummed coefficient function with the exact 0(a>s) 
one, in order to obtain a good approximation in the whole kinematical range. We can perform the 
matching in A^-space and then invert the final result to x-space; alternatively, we can invert the 
resummed coefficient function as in (14.261) . and then match it to the NLO x-space result. Given the 
high accuracy of our approach and the delicate issue of the exact Mellin transform, we have matched 
resummed and NLO coefficient function in both N- and x-spaces, and checked the consistency of 
our results. We discuss the matching in both spaces. 
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4.2.1 iV-space 



We would like to write the NNLL-resummed coefficient function as: 



Cres 
N 



+ 4 C) 



/ 2 \ 2 2 



rrir 



A 



(CO 

A' 



as(/4);/4>/4, m 



^s{^R),^ 2 R,^ 2 F ,m 2 z 



(4.27) 



where: 



1. K^°' [ols(hr)-> Hr, /j, F ,m z ] is a hard factor, introduced for the sake of including subleading 
terms, which corresponds to the difference between the constant terms which are present in 
the exact NLO coefficient function and the ones contained in the 0(as) expansion of the 
resummed result; 

2. Ajy [x; as(fi R ), n R , i^ F ,m z ] is the resummed coefficient function, presented in (14.51) : 

3. [as (/Mr), Hr, fJ/p, m z ] is a remainder function, collecting the left-over NLO contributions, 
which are suppressed at large N. 



The hard factor reads: 



K (C) 



<*s{iA),Vr,Vf,™ 2 z 



1 + 



«s(/4) c f ~„ 



7T 



with 



7T 



> + " 3, 
H F 3 



(4.28) 



(4.29) 



where Q(fj, F ,m z ) and Q\fi F) m z ) have been defined in Eqs. (14.31) and (14.251) . respectively. If we 
had used the step-function approximation, no constants would have been resummed, and the hard 
factor would have contained all the constant terms present in the NLO coefficient function, i.e. 
Q(/4,m|). 

(C) 

In iV-space, the remainder function d N is obtained subtracting from the exact NLO coefficient 
function in Mellin space [11] the 0(acs) expansion of K^A^ C \ We obtain: 



d 



(C) 
N 



as(/4)>/4>/4>"4 



/~i ( 2 \ 2 2 2 

[a s {fiR),fi R ,fi F , m z 
K ic) \a s (jJ? R ),ti%,t&>™l 



A 



(C) 
N 



as(/4),/4,/4, m ll) 

J / as 



R) 



c f 



71 



2^ (1) (iV) + ^(- + 



1 , mi 
— In— ^ , 

2 /4 \N 



+ 1 NJ T v \N + 1 jS 



If 



+ 



+ 



2 \N iV+1 
1 



+ 



7 



4iV 4(iV + 1) 



N 2 (N+iy 



(4.30) 
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4.2.2 x-space 



Likewise, we wish to write the NNLL-resummed coefficient function in x-space in a form analogous 
to Eq. (14271) : 



«s(/4)>/4>/4 



,m z 



K (c) [as(/^e),/4,/4> m ll A (c) far; a s (/4), fi% /4> m% 



+ d {c) x; a 5 (/4),/4,/4 



,m z 



(4.31) 



To get the remainder function and the constants, the factorized expression (12.131) of the coefficient 
function turns out to be particularly useful. First, we need the 0(as) expansion of the resummed 
result in x-space, that can be read from the integrand function in the Mellin transform (14.171) : 



(A {c) \x;as(ji%),ti%ti F ,m%]) = 5(1 - x) + 

\ L J /as 



«5(/i|)Cf 



7T 



1 , m^fl — x) 
In - z 



x 



4[1 - x] 



(4.32) 

Then, by comparing Eq. (I4.32[) with Eq. (12.131) . we obtain the constants K^ c \ i.e. the coefficient of 
the term ~ 0:5 5(1 — x) in Eq. (12.131) . obviously equal to the ones in Eq. (I4.28p . and the remainder 
function: 



x; a s (n 2 R ), n\,n 2 F ,m\ 



C [x;a s (/4),/4,/4,m| 
K (c) \a s (n 2 R ), 11% fi 2 F ,m z ] A {c) \x; a s ((i 2 R ), (i 2 R , fi F 



,m z 



a s (n 



7[ 



iJ C F {^(5 - -ir) - -(I -h.r) 



ln(l - x) + ln^f 

/4 



l+x\ \ 



x 



(4.33) 



In (I4.33j) . we denoted by C [x; as(fi 2 R ), fi R , H 2 F) m z ] the NLO MS coefficient function (12. 7L Of course, 
the Mellin transform of Eq. fT4T33l yields Eq. fl43U|) . 



5 Large- £ resummation in the initial condition 

In this section we resum in the NNLL approximation the threshold logarithms which appear in the 
initial condition of the perturbative fragmentation function, and we match the resummed expression 
with the exact NLO one. As discussed in [14], large-x resummation in the initial condition is 
pro cess-indep endent . 
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5.1 Resummed initial condition in Mellin space 



The initial condition (13. 4p present terms ~ [ln(l — x)/(l — x)} + and 1/(1 — x) + which are to be 
resummed to all orders. In iV-space, the Mellin transform of the NLO initial condition in Eq. (13.41) 
exhibits single and double logarithms of N: 



where the constant terms are given by: 



OFi m b 



l 0R) 
IT 



In 2 N + In 



rag 



2 lE + 1 In AT 



OFi 



1 



7l~ 

~6 



+ IE ~ 1e + \1E 



hi 



m b 
^OF 



The resummed initial condition has a generalized exponential structure [14], 



(D) 



M/4m)>/4m>/4^ m &] = exp{G^ } [a s (l4 R ), l4 R , Hof, m b] } 



(5.1; 



(5.2) 



(5.3) 



where 



G 



(D) 
N 



( 2 \ 2 2 2 
a S(HoR)> ^ORi ^OFi m b 



JV-1 



dz- 



1 



dk 2 



Z \Jm\{\-zf k 2 



■A 



a s {k 2 



D 



a s 



with k 2 and z defined as in (14. 5p . As in the coefficient function, the LLs in the exponent Gffl are 
~ agln n+1 N, the NLLs ~ a%\n n N, and so forth. To NNLL accuracy, we need A®, and A^\ 
given in the previous section, and the first two coefficients of 



n=l 



namely 



D (D 

D (2) 



-c F 

C F 



CA [l08-4 a3) + — + 54 



(5.5) 

(5.6) 
(5.7) 



Function D(as), called H(ots) i n [14], is characteristic of the fragmentation of heavy quarks and 
resums soft-gluon radiation which is not collinear enhanced. Its 0(as) coefficient can be found in 
[14], while can be read from the formulas in [27]. Moreover, as discussed in [38], D(ots) coincides 
with the function which resums large-angle soft radiation in heavy-flavour decays [23,39,40]. It is 
also equal to function S(as), which plays the same role in top-quark decay [16] and massive Deep 
Inelastic Scattering [41]. 
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The coefficient depends on the renormalization condition on the b mass and Eq. (15.7ft gives 
its value in the on-shell scheme. If rhj, is the 6-quark mass in another scheme, related to the pole 
mass nib via 

[i + 2£ fcWl 



7T 



(5i 



the following relation has to be fulfilled, with the coefficients A® clearly unchanged: 



dk 2 



m?(l-z) 2 k 2 



A 



as(e)} + D{as(ml(l-zf)} = A 

n Jm (l-z) 2 k z 



By solving the above equation order by order, one obtains the coefficients D®: 

f)W = D W. 
m = + 2k^A^. 



a s (k 2 )] + L)[a s (m 2 b {l - zf 

(5.9) 



(5.10) 
(5.11) 



As in the case of the coefficient function, we identify the integral ~ dk 2 A[as{k 2 )]/k 2 in Eq. (15.41) 
as a generalization of the NLO integral in Eq. ( 13.111) . Function D(as) generalizes the single loga- 
rithmic term, ~ 1/(1 — x) + , that we found in the NLO computation (I3.10p : 



- C, 



a s 



71 



D 



(i) ^s 



D 



71 



a s (mf(l - xf 



(5.12) 



Following the arguments discussed in Subsection 14.11 the Mellin transform in Eq. ( 15.41) will be again 



performed exactly and the term 



JV-l 



1) will not be approximated. The subleading terms 



which are resummed in this way can be obtained after expanding Eq. (15. 3p to O(as)'- 
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(D) 
N 



( 2 \ 2 2 2 



1 + 
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or) 
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N-l 



A (l) 
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(5.13) 
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Si(iV-l) -S 2 (iV- 1) 



- S 2 {N-1) 



and taking its large- iV limit: 

/ a (d) r / 2 \ 2 2 2 
[A N \a s {n 0R ),[i 0R ,n 0F ,m b 



1 + 



as 



71 

™2> 



+ Y\^mi) + O 



In 2 N + 1 - 2-i E - In 



N 



^OF 

ml 



(5.14) 



In AT 



(5.15) 



4 Since m& depends in general on a renormalization scale /j, m , one could estimate its contribution to the theoretical 
error on the cross section by varying — in additional to the usual factorization and renormalization scales — also 
/i m around m& inside a conventional range. 
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with 



OFi m b, 



mi I 



7T 



(5.16) 



5.2 Matching of resummed and NLO initial condition 

We follow a matching procedure analogous to the one for the coefficient function in order to obtain 
a reliable result throughout the full N (x) range. 



5.2.1 iV-space 



We would like to write the NNLL-resummed initial condition of the perturbative fragmentation 
function as: 



D 



N 



K (D) 



A 



(D) 
N 



The multiplying factor is obtained subtracting from the constants terms that are present in the 
NLO result, the ones which have been resummed: 



a s(t*l R ), Vim V<lF, m b 

(5.17) 



K {D) a s (nl R ),nl R ,nl F , 
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where 



''OF' m bJ 



Y"(fi 2 0F ,m 2 b ) = Y(f, 2 F ,m 2 b ) - Y'(fi 2 
and Y(fj,l F ,ml) and Y\n1 Fl m 2 ) given in Eqs. (15.21) and (15.161) . 
The remainder function for the initial condition reads: 
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(5.19) 
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(5.20) 
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.AT ' AT + 17 2N 2(JV+1) N 2 (N+l) 2 
Multiplying the NNLL-resummed coefficient function by the DGLAP evolution operator and the 



NNLL initial condition, we obtain the moments of the 6-quark energy distribution in e + e 



bb: 
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(5.21; 
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5.2.2 x-space 



The 0{ots) expansion of the resummed initial condition in x-space can be read from the integrand 
function of Eq. H5. 131) : 



(A {D) \x; a s {nl R ), fi 2 0R , fi 2 0F ,m 2 z ]) = 6(l-x)-\ 



<*s(j4r) c f 



7T 



In 



^OF 



1 — x rn 2 (l — x) 2 



(5.22) 

The constant is the coefficient of the ~ as 5(1 — x) term in (13.101) . obviously equal to (I5.18p . 
The remainder function can be as well obtained from Eq. (13.101) and reads: 



d {D) \x\ a s (jJ% R ), l4 R , li% F ,m\ 
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,m b 



as 

(5.23) 



On can check that the Mellin transform of Eq. (I5.23P agrees with Eq. (I5.20p . The NNLL- resummed 
initial condition in x-space is finally given by: 



D 1 



K (B) 



5 VoRi ^OF) m b 
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(5.24) 



6 Effective coupling constant 

In this section we introduce an effective QCD coupling constant (i) having no Landau pole and (ii) 
resumming absorptive effects related to parton branching to all orders. 



6.1 Space-like coupling constant 

If we denote by g M the typical 4-momentum entering in the renormalization conditions for the QCD 
coupling constant 0, the standard LO coupling constant reads: 

«5,lo(-<? 2 ) = ^[(-q 2 - it) / K 2 } = J Q [\n(\q 2 \/k 2 ) - mQtf)Y (0) 

There is a minus sign in front of the momentum squared q 2 , because of the opening of decay channels 
for the gluon (g —>■ qq, g —>■ gg) in the time-like region q 2 > 0. In order to have a renormalized real 



5 One can consider, for example, the symmetric point p 2 = p? = p 2 = q 2 in the qqg correlation function. 
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as, one usually considers a space-like configuration of the reference momenta: q 2 < 0. To avoid 
explicit minus signs, the LO expression of as is usually written as: 

as > LoiQ2) = /? ln(QVA 2 )' (6 ' 2) 

where Q 2 = — q 2 > in the space-like region. The specific properties of the QCD coupling 
constant involved in soft-gluon resummation in e + e~ annihilation are related to: 



1. an integration up to small momentum scales of the coupling constant, because of multiple 
parton radiation, giving rise to sub-jets with arbitrarily small masses. In the resummed 
expressions, (14.51) and (15.41) . one can indeed observe that the scale of as approaches zero once 
x — > 1. 

2. the kinematical configurations are always time-like, as we are considering multiple emissions 
in the final-state of e + e~ processes. 



Let us deal with the above issues by discussing first the analyticity properties of the standard 
coupling constant (I6.2p . as(Q 2 ) exhibits a cut for Q 2 < 0, associated with the branching of a 
time-like gluon (q 2 > 0) into 'physical' states, and the Landau pole for Q 2 = A 2 . While the 
former singularity has a clear meaning, the Landau pole is not physical and it just reflects the 
unreliability of the perturbative expansion for Q 2 pa A 2 . 

It was therefore suggested than one can replace the usual expression (16.21) with an analytic 
coupling as{Q 2 ), which has the same discontinuity as as(Q 2 ) along the cut Q 2 < 0, but is analytic 
elsewhere in the complex plane. As in [20-22], we write the analytic coupling constant as{Q 2 ) 
using the following dispersion relation: 

a s{Q 2 ) = 7T~ Disc s a g (-g), (6.3) 

2m Jo s + Q z 

where the discontinuity is defined as: 

Disc s F(s) = lim [F(s + ie) - F(s - is)] . (6.4) 

e — > 0+ 

Eq. (16.31) holds for Q 2 > 0, i.e. in the space-like region q 2 < 0: as in [21], we shall refer to it as our 
'space-like' analytic coupling constant. For Q 2 < 0, the integrand function in Eq. (16.31) presents a 
pole in the domain of integration. However, we can still give sense to Eq. (I6.3P for negative values 
of Q 2 , introducing a small imaginary part: Q 2 — > Q 2 + ie. 

Inserting in (16.31) the LO expression (16.21) . we get the LO space-like coupling constant: 

1 A 2 



«s(Q 2 ) = )r 



ln(Q 2 /A 2 ) Q 2 -A 2 



(6.5) 
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In (16.51) the Landau pole has been subtracted by a power-suppressed term, relevant at small Q 2 and 
negligible at large Q 2 ^> A 2 , where as(Q 2 ) still exhibits the same behaviour as as(Q 2 ). Likewise, 
including in the integrand function of (16. 3p the higher-order expressions of as(—s), one can get the 
space-like as(Q 2 ) to higher accuracy. 

6.2 Time-like coupling constant including absorptive effects 

Turning back to our calculation, we have resummed soft and/or collinear multiple radiation in the 
final state of e + e~ annihilation, i.e. a time-like parton cascade. Also, in our resummed expressions, 
(14.51) and (15.41) . the coupling constant is evaluated at a scale k 2 , which is roughly the transverse 
momentum of the emitted parton with respect to the radiating one. In Ref . [42] , it was in fact shown 
that, in the framework of resummed calculations, the momentum- independent coupling constant is 
to be replaced by the following integral over the discontinuity of the gluon propagator: 

as — > — / ds Disc, ^—^ — -. (6.6) 
2tt Jo s 

The integral (16. 6p is typically performed neglecting the imaginary part, ~ ire, in the denominator 
of as(—s) (see for example Eq. (16.11) ). i.e. assuming 

lnj^»7r (6.7) 

in the integrand function of ( 16.6[) . As a result, the integral ( 16. 61) turns out to be approximately 
equal to as evaluated at the upper integration limit: 

^- ( k ds Disc, ^tA ~ as (k 2 ). (6.8) 

27T JO S 

In our analysis, we wish to go beyond the assumption (16.71) and account for the terms ~ in in the 
denominator of as', this way, as pointed out in [23], one includes absorptive effects due to gluon 
branching, that are important especially in the infrared region. The new feature of our model is the 
fact that we avoid the Landau pole in the integral by using in the integrand function the space-like 
analytic coupling constant as(-s), just defined in Eq. (16. 3ft . 

As a result, our model consists in using the following 'time-like' effective coupling constant: 

a s (k 2 ) = [ k ds Disc, as< <~ S \ (6.9) 
2n Jo s 

In fact, Eq. (16. 9p makes sense only for k 2 > 0: the above integral would be zero for negative values 
of k 2 . We also remark a difference in our notation: the effective coupling ds(k 2 ) is function of the 
square of a four-momentum k 2 \ the standard as(— q 2 ) in Eq. (I6.2p is instead function of minus a 
squared four-momentum. At large k 2 , as(k 2 ) will be roughly equivalent to the standard as{k 2 ); 
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at small k 2 it will include non-perturbative power-suppressed effects. The goal of this paper is 
precisely to investigate whether including non-perturbative corrections using Eq. (16.91) everywhere 
in our calculation is suitable to reproduce the experimental data on £>-hadron production, without 
adding any further hadronization model. 

Inserting Eqs. (I6.5P in the integrand function of (16.91) . we obtain the LO time-like analytic 
coupling constant: 



as,ho{k 2 



2nip 



k 2 



k 2 



In In — + in] — In In — — in 



A 2 



A 2 



(6.10) 



Likewise, starting from the NLO «s(Q 2 ) in Eq. (I3.8p . we obtain the NLO time-like coupling con- 
stant: 



as,NLo(k 2 ) = a s ,Lo(k 2 ) + 



Pi 2ni 



In [ln(fc 2 /A 2 ) + in] + 1 In [\n(k 2 /A 2 ) - in] + 1 



ln(A; 2 /A 2 



At NNLO, the standard coupling constant reads: 

1 



in 



ln(P/A 2 



I7i 



.(6.11) 



a s (k 



] __ Pi ln[ln(fc 2 /A 2 )] /5 2 ln 2 [ln(A; 2 /A 2 )] - ln[ln(A: 2 /A 2 )] - 1 



+ 



/3 ln(P/A 2 ) 1 (3 2 ln(A; 2 /A 2 ) $ 

1 1 



ln 2 (fc 2 /A 2 ) 



/3 3 ln 2 (P/A 2 ) J ' 
and its time-like analytic counterpart: 



(6.12) 



&S,NNLo(k 2 ) 



+ 



as,m,o(k 2 ) - 
Pi ~ P0P2 \ 



P\ \ In 2 [ln^/A 2 ) + in] In 2 [ln(A; 2 /A 2 



m\ 



4ttz/3 5 [ [ln(F/A 2 



ITX\ 



[ln(A; 2 /A 2 ) - in] 



4niP% 



[ln(fc 2 /A 2 ) +tn] z [ln{k 2 /A 2 



in\ 



(6.13) 



In (16.131) we have also included the third coefficient of the /^-function: 



64vr 3 



2857 



C 



(^c 2 A + ^c A c F -cl ]nf + 



79^ 11 , 

— C A H C F n 2 

54 9 



(6.14) 



In Eqs. (16.101) . (16. lip and (16 . 1 3[) we have used a complex notation, which yields quite compact 
expressions for the time- like coupling constant. However, we can rearrange the above equations 
and express as{k 2 ) as a real function of k 2 , as done in [39]. 

Before closing this section, we point out a few more issues. Expanding ots{k 2 ) for ln(/c 2 /A 2 ) 3> n, 
it is possible to relate standard and analytic time-like coupling constants: 



a s {k 2 ) = a s {k 2 ) - 



{kPo) 



aUk 2 ) + 0{a%). 



(6.15) 
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We also need to modify the matching condition for the strong coupling constant when running from 
rif to Uf-x active flavours. In terms of the standard as(k 2 ), it reads [43]: 

o%)K) = as(n r i)(mp - T^<4(n f -i)(™i) + °( a s), ( 6 - 16 ) 

where fh q is the running quark mass in the MS renormalization scheme, i.e. fh q = m^ s (m q ). When 
using the analytic effective coupling a s , we shall have to modify Eq. (16.161) according to Eq. (I6.15p . 
observing that (3q, given in Eq. (13.71) depends on rif. We obtain: 

^S(n } )(m 2 q ) = as {ni -i){rn 2 q ) - ~ ^^^ j "1(^-1) + °( a s)- ( 6 - 17 ) 

Of course, Eqs. (16.161) and (16.171) can be also expressed in terms of the pole quark masses. 



7 Modelling non-perturbative corrections 

In this section we describe our model for non-perturbative effects in bottom-quark fragmentation, 
based on the effective QCD coupling constant considered above. Some properties of the model have 
already been anticipated in the previous section; here we present a systematic discussion. 

We shall use the effective coupling constant as(k 2 ) defined through Eq. (16.91) in place of the 
standard one, in order to include power corrections. The dominant non-perturbative effects in 
6-fragmentation occur for 

1- — <ar<l, (7.1) 
m b 

where \i is of the order of the QCD scale, i.e. /i ~ O(A). Within our model, such contributions 
are associated with soft interactions of the b quark in the fragmentation into a B hadron and are 
analogous to the well-known Fermi motion of a decaying b quark inside a B. In the perturbative 
fragmentation approach, such effects become relevant when we evaluate at large x the quantities 

C = m z \/\ — x and S = rrib(l — x), (7.2) 

whose squares are the limits of the dk 2 /k 2 integration in the resummed coefficient function (14.51) 
and initial condition (15. 4p . respectively, as well as the scales of functions B(as) and D(as)- It is 
interesting to evaluate C, S and the corresponding values of as for x = 0.8, since, as we shall show 
in the following section, the 5-hadron spectrum in e + e~ annihilation is peaked about this value of 
x. We findS 

C ~ 40 GeV , a s (C 2 ) ~ 0.13 ; 5-1 GeV , a s (S 2 ) ~ 0.33. (7.3) 
6 For x = 0.9, the corresponding numbers are: C ~ 30 GeV, a s (C 2 ) ~ 0.14, S ~ 0.5 GeV, a s (S 2 ) ~ 0.44. 
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It is therefore clear that non-perturbative effects are more relevant in the initial condition (15. 4R . 
where S plays a role, rather than in the coefficient function (14.5j) . depending on C. In order to deal 
with such effects, we shall insert as in place of the standard as in the resummed initial condition 
(15.41) and, as part of our model, we shall perform the Mellin and inverse Mellin transforms exactly. 
As discussed in Subsection 14.11 the issue of the integration over z when using the effective coupling 
is delicate and does deserve a deeper investigation in the next future. However, it was pointed 
out in [44] that, in order to include the power-suppressed corrections originated by the effective 
coupling, ~ 0(1 /N) or ~ 0(1 — x), performing the ^-integration exactly is necessary. In fact, an 
approximate Mellin transform, using the step function (14. 16[) or the formulas beyond NLL in [33], 
would suppress most of such effects [44]. The Mellin transforms were performed exactly even in 
Ref. [9], where the effective coupling was used to describe non-perturbative effects in 5-meson 
decays. 

In the coefficient function, the use of as(k 2 ) and the exactness of the Mellin transform are less 
crucial than in the initial condition of the perturbative fragmentation function. In any case, for 
practical convenience, we shall still use the effective coupling constant and perform the z-integration 
exactly. Besides, as(k 2 ) will be employed everywhere in our calculation, including the constant 
terms and the remainder functions presented in Eqs. (14.301) . (14.331) . (I5.20p and (I5.23p . 

We also remark that that there is not a unique way to construct a model based on the analytic 
coupling constant. In fact, we have to make two choices: 



1. we can include the absorptive effects, which are always present in time-like kinematics, and 
use as{k 2 )] alternatively, we do not include such effects and employ the space-like ag(fc 2 ); 

2. we can perform a power expansion of the higher orders or not. 

The latest point deserves some more comments. By 'power expansion' we mean that higher 
orders have in front a power of the effective coupling constant: 



a n s (k 2 



— / ds Disc, 
2% Jo 



a s [-s 



(7.4) 



The n-th power is taken after the discontinuity of as{—s) is computed and the integral over 
the gluon virtuality s is performed. Adopting the 'non-power expansion' choice, as originally 
proposed in [21], consists instead in evaluating first the discontinuity of ot s (— s), and then 
integrating over s. Formally, ag(fc 2 ) will have to be replaced according to: 



a n s {k 2 ) -> a^k 2 ) = 



2tt Jo 



ds Disc, 



(7.5) 



Unlike ( 17.41) . Eq. ( 17. 5p is linear in the discontinuity. 



We have followed an empirical criterion, leaving to the future the task of a theoretical justifi- 
cation: we select the model which gives spectra closer to the data, without adding any further 
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non-perturbative fragmentation function. We can anticipate that have found that the power ex- 
pansion of the coupling constant which includes the absorptive effects, i.e. as(k 2 ), leads to the best 
description of the experimental data. To be more general, we stress that our choices in modelling 
non-perturbative effects are partly related to the accuracy of the perturbative resummed calcula- 
tion that we are using. Considering, for example, function v4[a,s(fc 2 )], we have used its expansion 
to the third order and the NNLO as{k 2 ). However, if we were to know function A to a different 
level of approximation or, in principle, even to any order, a different effective coupling constant 
may still yield the same A[as(/c 2 )] and the same 5-hadron spectrum. 

Since Eq. (16. 15ft is basically a change of scheme for the coupling constant, the coefficients of 
the terms of 0{a%) in functions A(ats), B(as) and D(as), appearing in the resummed expressions 
(14.51) and (15.41) . are to be modified. In our NNLL approximation, we need to replace A^ according 
to: 

A® -> A® = A ® + M>V>, (7.6) 

3 

where we have made use of Eq. (I6.15p . The other coefficients entering in functions A(cts), B(a>s) 
and D(ocs) at NNLL are instead left unchanged when replacing as{k 2 ) with as{k 2 ). 

Let us now discuss another delicate point of our model. In standard resummations, which use 
the standard coupling constant and resum only the logarithms of the Mellin variable N, there is a 
clear counting of the terms in as{k 2 ) which are to be kept or dropped. For example, in the NLL 
approximation: 



.4 

That is because 



a s (k 2 )} -> A« as ' NLo(fc2) + A^ 



TT TT 2 



(7.7) 



a ^ {k2) ~ Tnwwy (7 - 8) 

and therefore it has, for k 2 — > oo, the same asymptotic behavior as cts,NLo{k 2 ), which contains a cor- 
rection term behaving like ln(ln(fc 2 /A 2 ))/ In (k 2 /A 2 ) ~l/ln (k 2 /A 2 ). In the standard resummation 
scheme — which is a kind of minimal scheme — only logarithms are resummed and the asymptotic 
expansion of the coupling constant is used to organize the series of the infrared logarithms. On 
the other hand, the large-fc 2 expansion of our time-like as{k 2 ) is much more complicated; even the 
lowest order as,ho(k 2 ), proportional to l/Po, presents terms of any order for k 2 — ► oo: 

as ' Lo{k2) ~ + lu 3 (P/A 2 ) + \n\k 2 /A 2 ) + " " (7 ' 9) 

In fact, we observe that: (i) in our model we would like to include as many contributions as possible; 
(ii) there is no real reason to neglect higher-order corrections to as(k 2 ), proportional to Pi, (3 2 an d 
so on. Therefore, we find it safe using «5 i nnlo(^ 2 ) everywhere in the resummed expressions. For 
example: 

~ . A^ A^ _ i( 3 ) ^ 

A[d s (k 2 )} -> — a s ,NNLo(k 2 ) + — ^ a 2 SNNLO (k 2 ) + — T a 3 SNNLO (k 2 ) . (7.10) 
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In this way we include all the logarithmic terms of the standard NNLL resummation, plus some 
subleading contributions. Another possibility might have been: 

4(1) 4(2) ^(3) 

A{d s {k 2 )) -> a5,NNLo(fc 2 ) + — a|,NLo(fc 2 ) + — «I,lo(^)- (7-11) 

However, at small k 2 , where we are mostly sensitive to non-perturbative effects, higher-order correc- 
tions to the time-like coupling constant are negative and sizable. Comparing, in fact, Eqs. (I6.10p . 
(16. lip and (16.131) for low values of k 2 , we find: 

«S,NNLo(^ 2 ) < aS,NLo(^ 2 ) < ®S,Lo(k 2 ). (7.12) 

Hence, if we used Eq. ( 17.1 ip rather than Eq. (I7.10p . the term proportional to would be enhanced 
at large x. We can anticipate that this would worsen the comparison with the experimental data. 



8 Phenomenology — x-space 

In this section we present results on the B-hadron energy spectrum using the resummed partonic 
calculation based on the perturbative fragmentation formalism, and modelling non-perturbative 
effects by means of the time-like effective coupling constant. The basic assumption of our model 
is that, whenever we use ots{k 2 ) instead of as{k 2 ), the 6-quark energy fraction x will have to be 
replaced by its hadron-level counterpart xb- 

2pB ' q (8.1) 



x B 



with pb being the momentum of a 6-flavoured hadron produced in e + e annihilation. 



The .B-hadron spectrum in moment space can be obtained from Eq. ( 15.211) . after replacing the 
standard coupling constant with the analytic time-like one: 
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(8.2) 



In Eq. (18.21) we assume that the coefficient function and the initial condition are resummed in the 
NNLL approximation, and that the NNLO analytic coupling constant (I6.13P is used. In order to 
recover the x-space results, we perform the inverse Mellin transform of Eq. (18. 2p : 

'- -N ^( B )^ 2 ,,2 ,,2 ,,2 ,,,2 ,,,2 



(B) ( 2 2 2 2 2 2 

a \Xb, Hr, Hor, flQFy f^Fi m bi m Z 
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^ 2 2 2 2 2 \ 
L Ri ^ORi ^OFi A*F) m bi m z)i 



•3) 



where c is a positive constant. We point out that the inverse transform (18.31) is computed in the 
standard mathematical way: no prescription to deal with the Landau pole, such as the minimal 
prescription [36], is needed. That is because our effective coupling as(k 2 ) does not present the 
Landau pole any longer. Since the inverse transform is made in a numerical way, we have checked 
that the results are stable with respect to change of the integration contour, i.e. of c. 
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8.1 5-hadron spectrum and comparison with experimental data 



Following [18], we consider data from SLD [1], ALEPH [2] and OPAL [3] collaborations at the 
Z° pole. ALEPH reconstructed onlvi? mesons, while the SLD and OPAL samples also contain a 
small fraction of 6-flavoured baryonsLj. In principle, one should consider such data separately, as the 
hadron content is different; also, from the theoretical viewpoint, as pointed out in the introduction, 
there is no real reason why the same model should describe both meson and baryon production. 
However, the SLD and OPAL samples are inclusive and the baryons are anyway very little. In 
fact, Ref. [18] found that it is possible to describe all data points fitting the Kartvelishvili non- 
perturbative fragmentation model [5] or the cluster and string models implemented in HERWIG [45] 
and PYTHIA [46]. It is therefore reasonable using the same hadronization model, in our case the 
effective coupling constant, for the comparison with all three experiments, and investigating how 
it fares against all data and against each experimental sample. As in the previous analyses, when 
doing the comparison, we neglect the correlations between the data points and sum the experimental 
systematic and statistical errors in quadrature. 

Unlike Refs. [15-18], we shall not perform a fit to the data, since we do not have any tunable 
parameter in our model, apart of course from the ones contained in our perturbative calculation. 
Rather, we shall investigate the theoretical uncertainty on our prediction, by varying the parameters 
in our computation, such as scales and quark masses. The default values of our parameters will be: 

= = rn z ; fi m = /i 0F = m b , (8.4) 

where and //p are the renormalization and factorization scales in the coefficient function (14.271) . 
respectively, and fion and hof in the initial condition (15.17 ). Consistently with the Particle Data 



Group [47], we set mz = 91.19 GeV and as(m 2 z ) = 0.119^. For the MS quark masses, entering in 
the matching condition (I6.16p . we choose: rh b = 4.2 GeV, rh c = 1.25 GeV, rh s = 0.1 GeV, while 
the up- and down-quark masses will be neglected. The choice of rrif,, the 6-quark mass entering 
in the initial condition of the perturbative fragmentation function, deserves some extra comment. 
In fact, using the pole or the MS b mass in the initial condition is irrelevant at NLO and when 
resumming threshold logarithms up to NLL accuracy. Beyond such approximations, the choice of 
the renormalization scheme makes a difference. Refs. [27,28], which calculate the initial condition 
to NNLO, use the pole mass in their computation. Although in the matching we are using the NLO 
initial condition, we are still relying on Ref. [27] for the coefficient in the NNLL resummation 
(15.41) . and therefore we should use the pole mass as well. However, since we are aiming at predicting 
the I?-hadron energy distribution, it is not uniquely determined whether, after using the analytic 
coupling constant to include power-suppressed effects, m b should be the 6-quark or the S-hadron 
mass. We believe that it is safe adopting a quite conservative choice, i.e. 4.7 GeV < m h < 5.3 GeV, 
which includes the present estimations for the b pole mass as well as 6-flavoured hadron masses [47]. 
Our default value will be nib = 5 GeV. 



7 A naive estimation based on the 1/N C expansion, where N c is the colour number, would predict a fraction of 
about l/iVf « 10% 6-flavored hyperions compared to B mesons. 

8 This choice corresponds to as{m 2 z ) = 0.117, as can be obtained from Eq. (|6.15|) . 
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Figure 1: U-hadron spectrum yielded by our resummed calculation using the analytic time- like coupling 
constant as(k 2 ) to model non-perturbative effects. We investigate the dependence on the factorization 
scales. Solid lines: //of = Tn b /2, m b and 2m;,; dashed lines: (ip = fnz and 2mz- The other 

parameters are kept to their default values. 

Before presenting our results, we point out that this kind of systematic analysis of the theoretical 
uncertainty was not performed, e.g., in Refs. [16-18]. In fact, when using a hadronization model, 
the fitting procedure would possibly adjust the free parameters to reproduce the data, even when 
varying the inputs of the parton-level computation. As long as a reliable hadronization model 
is used, changing the perturbative quantities will just lead to different best-fit parameters in the 
non-perturbative fragmentation function. 

In Fig. [Tj we show the data points and our prediction, investigating its dependence on the 
factorization scales \Lp and //of- in Fig.[2]we look instead at the dependence on the renormalization 
scales fiR and [Iqr. For //# and //p we choose the values mz/2, mz and 2m^; for // h and //of, 
m b /2, nib and 2m b . We vary each scale separately, keeping all other quantities to their default 
values, in order to avoid an excessive number of runs. In Fig. [TJ we notice that the dependence of 
our prediction on fi p, the factorization scale entering in the initial condition of the perturbative 
fragmentation function, is rather large, especially at small xb and around the peak. At small xb, 
the spectrum varies up to a factor of 2 if /j, Q p changes from m b /2 to 2m;,; around the peak the 
impact of the choice of //of is about 20%. The effect of the value chosen for \ip is instead pretty 
small, and well within the band associated with the variation of //of- A possible explanation of 
the fairly large dependence on //of could be the fact that, although we are working in the NNLL 
approximation, we are still matching the resummation to the NLO exact result, and not to the 
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Figure 2: Dependence on the renormalization scales fiR and /iqr. Solid lines: fj,R varied between 
and 2mz; dashed lines: fion between and 2m&. 

4 I — i — i — i — I — i 1 — i — i — I — i — i 1 — i — I — i — i — i 1 — I — i — i 1 — i — I — r 




Figure 3: Dependence of the U-energy distribution on as(rn%) and on rut,- Solid lines: standard coupling 
a s( m z) varied from 0.117 to 0.121; Dashed lines: b pole mass varied from 4.7 to 5.3 GeV. 



31 



NNLO one. This generates a mismatch between the NNLL terms in the resummed expressions 
(~ OghiN, a| In 2 iV, . . .) and the remainder functions, given in Eqs. (14.301) and (15.201) . which 
have been included only to NLO. This mismatch is more evident in the initial condition, where 
/i F plays a role, since the coupling constant is larger, being evaluated at scales that are smaller 
than in the coefficient function. We should expect a milder dependence on such a factorization 
scale if we matched the resummed initial condition to the exact NNLO result [27]. Moreover, the 
scale /io_F also enters in the DGLAP evolution operator (13.61) . which we have implemented to NLL 
accuracy, using NLO splitting functions. Ref. [26] has recently calculated the NNLO corrections to 
the time-like splitting functions, which include a contribution analogous to the one ~ entering; 
in the resummed expressions (14. 5p and (15.41) . The fact that we have not included such effects 
may be a further source of mismatch, which may be fixed implementing the splitting functions 
to NNLO accuracy. In any case, a peculiar feature of our model is that, since we are not using 
any hadronization model, our theoretical error includes, at the same time, uncertainties of both 
perturbative and non-perturbative nature. The partonic calculations in Refs. [14,16,17], which 
are NLL/NLO and use the standard coupling constant, yield indeed a very mild dependence on 
the quantities which enter in the perturbative calculation. However, in order to predict hadron- 
level observables, these perturbative computations need to be convoluted with a non-perturbative 
fragmentation function, whose parameters, after being fitted to SLD and LEP data, exhibit errors 
which are typically of the order of 10%, in such a way that the hadron spectra still present fairly 
large uncertainties. 

Turning back to Fig. [U let us observe that our distributions become negative and oscillate at 
very small and at very large xb- In fact, the coefficient function (12.71) presents a term of the form 

aM) ° F lux, (8.5) 

7T 

which is enhanced at small x and has not been resummed in our analysis. In any case, that does 
not affect much the comparison with the data, since the latter are at x > 0.12. The coefficient 
function (12.71) and the initial condition (13. 4p also contain terms of the type 

- ° F HI - x) and 2 as{ ^ )CF Hl-x), (8.6) 

7T 7T 

which become large for x — > 1 and have not been resummed. It is at present not known how to 
accomplish this task. The contributions in Eq. (18. 6p are responsible of the oscillating behaviour 
at large xb- It is conceivable that the inclusion of NNLO corrections to the coefficient function 
and the initial condition may partly stabilize the distribution at the endpoints. In any case, we 
are aware that our simple model, based on an extrapolation of the perturbative behaviour to small 
energy scales, is expected to fail for very large x. It is therefore safe discarding a few data points at 
large xb when comparing with the data and, e.g., limiting ourselves to xb < 0.92. We find that our 
default parameters give a good description of the ALEPH data (x 2 /dof = 21.4/16), but reproduce 
rather badly the OPAL (x 2 /dof = 162.7/18) and SLD (x 2 /dof = 109.1/20) ones. The overall 
X 2 /dof, computed as if all measurements were coming from one experiment is also quite large: 
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Table 1: Results of the comparison of our model with LEP and SLD data, for the most significant values 
of fiQF and m&. We have set oes(m z ) = 0.119, while the other parameters have very little impact. 





m h 


X 2 /dof (ALEPH) 


X 2 /dof (OPAL) 


rVdof (SLD) 


X 2 /dof (overall) 


m b 


5 GeV 


21.4/16 


162.18/18 


109.1/20 


293.2/54 


m b /2 


5 GeV 


24.8/16 


30.4/18 


47.8/20 


103.0/54 


m h 


5.3 GeV 


11.9/16 


116.1/18 


84.2/20 


212.2/54 



X 2 /dof = 293.2/54. A much better description of SLD, OPAL and the overall sample is obtained 
setting hqf = rn b /2. We obtain: x 2 /dof = 24.8/16 (ALEPH), 30.4/18 (OPAL), 47.8/20 (SLD) 
and 103.0/54 (overall). Such values of x 2 /dof are perfectly acceptable, since we are comparing 
data from different experiments and, as already discussed, SLD and OPAL, unlike ALEPH, also 
reconstructed a small fraction of 6-flavoured baryons. 

The theoretical uncertainties due to the values chosen for as(m z ) and mj are explored in Fig. [3j 
we consider as(m z ) = 0.117, 0.119 and 0.1210, and = 4.7, 5.0 and 5.3 GeV. The impact of 
the choice of these quantities is comparable and well visible throughout all the a^-range, though 
smaller than the one due to the variation of fi p. The effect is about 10% at average values of 
xb and grows up to 35% at large x. While changing as(m 2 z ) and nib does not improve much the 
comparison with OPAL and SLD, an excellent description of the ALEPH data (x 2 /dof = 11.9/16) 
is obtained for nib = 5.3 GeV, a value compatible with 5-meson masses, and still as(m 2 z ) = 0.119. 
Nevertheless, the comparison with the other experiments gets worse for this value of the bottom- 
quark mass, as we get x 2 /dof = 116.1/18 for OPAL, 84.2/20 for SLD and 212.2/54 overall. As for 
the MS masses nib, f^c and ni s , the effect of their variation on the x^-spectrum is very little: we 
do not show such plots for the sake of brevity. The most relevant values of x 2 /dof when we vary 
the perturbative parameters are collected in Table [TJ 

8.2 Relevance of NNLL effects 

Before turning to the analysis in moment space, we wish to estimate the impact on our prediction 
of the NNLL threshold-resummation corrections. 

In Fig. H] we present the experimental points and the spectra yielded by our model using the 
parameters which give the overall best fit to the data (/iqf = nib/ 2 and nib = 5 GeV), but 
within two approximations schemes: NNLL soft resummation along with NNLO effective coupling 
constant (solid line), and NLL soft resummation with NLO coupling constant (dashed). In the 
NLO/NLL prediction, we still perform the Mellin transforms of the resummed expressions exactly, 
i.e. without the step-function approximation (14.161) . By dropping the NNLL terms, i.e. the 
contributions proportional to and in the resummed expressions ( 14.5}) and ( 15. 41) . we 

9 The corresponding range for the analytic time-like coupling constant is 0.115 < ots(m%) < 0.119. 
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see that the spectrum gets shifted to larger xb, and the agreement with the experimental data 
gets worse 0- For x B < 0.92 we obtain: x 2 /dof = 11.18 (SLD), 7.01 (ALEPH), 16.27 (OPAL) 
and 12.05 (overall). The inclusion of the NNLL threshold contributions to the coefficient function 
and to the initial condition of the perturbative fragmentation function, together with the NNLO 
corrections to cts{k 2 ), is therefore crucial to get reasonable agreement with the data and the x 2 /dof 
quoted in Table [H In fact, when using the time-like coupling constant, the 0(ot\) coefficient 
in the NNLL-resummed coefficient function and in the initial condition gets enhanced according to 
Eq. ( 17. 6p . The change — > A^ shifts the peak of the spectrum towards a lower value of xb- 
If we had used the space-like coupling constant as{k 2 ), given in Eq. (16.31) . instead of as(k 2 ), we 
would have obtained the same A^ and a rather poor description of the data. 

It is also interesting to gauge the overall impact of the inclusion of power corrections by means of 
our model. In Fig. H] we also show the NLL-resummed parton-level spectrum of Ref. [14], where the 
standard NLO coupling constant is used, the Mellin transforms are performed by the step-function 
approximation, and the inversion to x-space is done using the minimal prescription [36] to avoid 
the Landau pole. We see that the distribution of Ref. [14] is peaked at very large x and is very 
far from the data for xb > 0.4. In fact, the 6-quark spectrum of [14] is expected to be convoluted 
with a non-perturbative fragmentation function, whose parameters need to be fitted to reproduce 
the data. Our spectrum is also broader than the one of [14] : one can actually show that this is due 
to the fact that we have performed the Mellin transforms exactly. 



9 Phenomenology — TV-space 

We would like to present the results yielded by our approach in Mellin moment space, and compare 
them with the measurements of the DELPHI experiment [4]. It was argued in [19] that working 
in iV-space is theoretically preferable, since one does not need to assume any functional form for 
the non-perturbative fragmentation function and can fit directly its moments. Moreover, as shown 
in [7,18], the moments obtained fitting a non-perturbative model discarding data points at small 
and large x B are quite uncorrect, since the tails of the distributions play a crucial role to obtain 
the right iV-space results. Our case is clearly different, since we are not tuning a non-perturbative 
fragmentation function, and therefore we do not have the problems related to the fits. However, 
our distributions are still negative and unreliable at very small and large xb, and we were forced 
to discard few data points to obtain acceptable values of % 2 /dof. It is therefore still interesting 
to present the moments obtained using the analytic coupling constant and explore the theoretical 
uncertainties, by changing the parameters which enter in the calculation, along the lines of our 
x-space analysis. We calculate the B spectrum in Mellin space directly from Eq. (18.21) . and choose 
the same default values for scales, quark masses and as(m|) as in x-space. 

10 Similar conclusions are reached by comparing thrust, heavy-jet mass and C-parameter distributions — computed 
with the effective timedike coupling within NLL accuracy — with LEPf data [48]. 
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Figure 4: The solid line is our best prediction for £?-hadron production (j = B), obtained using NNLL 
threshold resummation and the NNLO time-like coupling, super-imposed to the data. The dotted curve 
is obtained by dropping the NNLL terms, i.e. working within NLL accuracy. The dashed line is the 
parton-level spectrum of Ref. [14] (j = b), which uses the standard as(k 2 ) and NLL large- x resummation. 



In Table |2] we present the results of our study: we quote the experimental moments measured 
by DELPHI, the central values yielded by our calculation, denoted by (er^) t h, and the errors due 
to all the quantities which we vary. The data correspond to N = 2, 3, 4 and 5; the first moment 
is cr|Li = 1, since both data and our spectra are normalized to unity. We estimate the overall 
theoretical error summing in quadrature the errors due to the variation of each parameter. 

We observe that the experimental moments exhibit very little errors and that our central values 
are smaller by about 5 -j- 10% than the DELPHI ones. This result confirms what we had found 
in xs-space: our x^-spectra go to zero more rapidly than the experimental data for large xb, arid 
become negative for i B > 0.96; for low and intermediate Xb, our curves lie above the data. It is 
therefore reasonable that our moments are smaller than the measured ones. However, we find that, 
within the uncertainties due to masses and scales, our calculation is in fair agreement with the 
data. As observed for the a^-spectrum, the uncertainty due to the choice of fiQp is pretty large, 
and the ones due to as(m%) and rnf, are smaller but visible. The other scales and the MS masses, 
entering in the matching of as at different flavour numbers, have very little impact on the moments 
of the B cross section. As far as the uncertainties are concerned, we observe that the DELPHI 
data are for N < 5, hence the terms ~ ln fc N in the resummed exponents are not so dominant 
with respect to other contributions, such as the constants. Therefore, even in iV-space, we expect 
that the further inclusion of NNLO contributions will have a significant impact on the prediction 
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Table 2: Moments a§ from DELPHI [4] and moments [cr^] t h yielded by our calculation. We quote the 
uncertainties due to the parameters which enter in the perturbative calculations, varied as discussed in 
the text. The theoretical total error is estimated as the sum in quadrature of the partial errors. We also 
present the moments a b N according to Ref. [14]. 





(x) 


(x>) 


(x 3 ) 


(x 4 ) 


e + e data cr^ 


0.7153 ±0.0052 


0.5401 ± 0.0064 


0.4236 ±0.0065 


0.3406 ± 0.0064 




0.6867 ±0.0403 


0.5019 ±0.0472 


0.3815 ±0.0465 


0.2976 ±0.0462 




0.0014 


0.0011 


0.0009 


0.0007 




0.0066 


0.0067 


0.0059 


0.0051 




0.0022 


0.0028 


0.0031 


0.0033 


<5cr^ Oof) 


0.0364 


0.0414 


0.0398 


0.0364 


5a%(m b ) 


0.0111 


0.0145 


0.0153 


0.0150 


5a%(rh b ) 


0.0004 


0.0005 


0.0006 


0.0006 


<5cr^(m c ) 


0.0003 


0.0005 


0.0006 


0.0006 


5a§(rh s ) 


0.0004 


0.0007 


0.0008 


0.0008 


5a§(a s (m 2 z )) 


0.0113 


0.0158 


0.0173 


0.0176 


< 


0.7734 ± 0.0232 


0.6333 ± 0.0311 


0.5354 ± 0.0345 


0.4617 ± 0.0346 



and decrease the theoretical error. 

We also present in Table[2]the moments yielded by the NLL-resummed calculation of [14], which 
uses the NLO standard as(k 2 ), whose x-space results have already been displayed in Fig. HI We 
estimate the error on the moments of Ref. [14] as we did for the one on our calculation, varying the 
perturbative parameters in the same range. We find that the iV-space results of [14], denoted by 
a b N to stress that they are a parton-level result, are very far from the data and much larger than 
the moments obtained using the analytic coupling constant. This result is in agreement with the 
plots in Fig. HI and confirms the remarkable impact of non-perturbative effects in iV-space as well. 
Moreover, we observe that the errors on the moments of [14] are smaller than the ones yielded by 
our calculation. In fact, Ref. [14] employs NLL threshold resummation and DGLAP evolution, and 
matches the resummation to the NLO results. As already discussed in the x-space analysis, using 
NNLL large-x resummation but still NLL DGLAP evolution and NLO remainder functions, as we 
did, may lead to a mismatch which produces larger uncertainties. We expect that the inclusion 
of NNLO corrections to the initial condition and to the splitting functions will lead to a weaker 
dependence on the input parameters of our analysis. In any case, the perturbative moments of [14] 
need to be multiplied by the moments of a non-perturbative fragmentation function extracted from 
the data, which will lead to larger uncertainties on the hadron-level moments. 
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10 Conclusions 



We studied I?-hadron energy distributions in e + e annihilation at the Z° peak by means of a model 
having as the source of non-perturbative corrections an effective QCD coupling constant. 

The physical idea behind our model is that the main non-perturbative effects, related to soft 
interactions in hadron bound states, are not very strong coupling phenomena, but can be described 
by an effective coupling constant of intermediate strength, typically as ~ 0.3 — 0.5. That implies 
that, within our model, these non-perturbative corrections — the Fermi motion in B decays being 
the best-known example — can be obtained from perturbation theory by means of an extrapolation. 

We modelled power corrections using a NNLO effective a s {k 2 ) constructed removing first the 
Landau pole, according to an analyticity requirement, and then resumming to all orders the ab- 
sorptive effects related to time-like gluon branching. The resulting as(k 2 ) significantly differs from 
the standard as{k 2 ) at small scales. We described 6-quark production within the framework of 
perturbative fragmentation functions, resumming threshold logarithms in the coefficient function 
and in the initial condition of the perturbative fragmentation function to NNLL accuracy, and 
matching the resummed results to the NLO ones. The perturbative fragmentation function was 
evolved using NLL DGLAP equations. When implementing threshold resummation in iV-space, as 
part of our model, we decided to perform the Mellin transforms exactly, which leads to the further 
inclusion in the Sudakov exponent of terms ~ ln fc N of higher order (N 3 LL, N 4 LL, . . .), along with 
some constants and power-suppressed terms O (l/N). 

We presented results on the 5-hadron energy distribution, relying on this model to include 
power corrections, without introducing any non-perturbative fragmentation function with tunable 
parameters. When studying the B spectra, we investigated the dependence of our prediction on 
the quantities which enter in the perturbative calculation, such as quark masses, renormalization 
and factorization scales, and as(m 2 z ). We observed that the x^-distributions become negative and 
oscillating at very small and large xb, which is due terms ~ lnx and ~ ln(l — x), that are present in 
the NLO remainder functions and have not been resummed. We have therefore discarded a few data 
points at large xb, where our calculation is anyway unreliable, and limited ourselves to x B < 0.92. 
We compared our results with OPAL, ALEPH and SLD data on 6-flavoured hadron production 
and found that, within our theoretical uncertainties, our calculation is able to reproduce quite well 
the ALEPH and OPAL data, while it is marginally consistent with SLD. As for the dependence 
on the perturbative parameters, the effect of the choice of the factorization scale /iof, appearing 
in the initial condition of the perturbative fragmentation function, is quite large and the data are 
better described for low values of /iof- The dependence on asijn 2 z ) and on nib is also quite visible: 
for example, a value of nib consistent with S-meson masses, which is reasonable since our model 
assumes Eb — Eb, yields an excellent description of the ALEPH data. The other parameters have 
instead very little impact on the a^-spectrum. 

We also compared our results in Mellin space with the moments measured by the DELPHI 
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collaboration. We found that the central values of our iV-space results are smaller than the experi- 
mental ones, but nonetheless, within the theoretical uncertainties, our moments are consistent with 
the data. In particular, a fairly large uncertainty is still due to the choice of fj, p, as was observed 
in the rr-space analysis. 

In summary, we find it remarkable that, within our theoretical uncertainties, we have been able 
to describe the LEP data in both x- and iV-spaces, by modelling non-perturbative corrections via 
the analytic time-like coupling constant and without tuning any parameter to such data. It is 
also pretty interesting that a model which succeeded in reproducing the photon- and hadron-mass 
energy distribution in £?-meson decays [9] has led to a reasonable fit of £> -pro duct ion data, although 
they are quite different processes, characterized by different energy scales. 

Our model looks therefore quite promising and we plan to further apply it to other processes 
and observables. A straightforward extension of the study here presented is the investigation of B 
production in top and Higgs decays, using the perturbative calculations in [15-17], and modelling 
power corrections as in the present paper. We can also use our model along with Monte Carlo 
event generators: in Ref. [18], the hadronization models of HERWIG and PYTHIA were tuned to 
the same data as the ones here considered, and then used to predict B production in top and Higgs 
decays. We may thus think of using the parton shower algorithms of HERWIG and PYTHIA to 
describe perturbative 6-production, with our analytic coupling constant in place of the cluster and 
string models. At the end of the cascade, we can assume Ef, — > Eb and compare the results with 
the data. The spectra yielded by Monte Carlo event generators are positive definite and do not 
exhibit the problem of becoming negative. However, parton shower algorithms are equivalent to a 
LL/LO resummation, with the further inclusion of some NLLs [49]. Since we have learned from 
this analysis that the inclusion of corrections of higher orders is crucial to reproduce the data, we 
may have to add to the Monte Carlo Sudakov form factor contributions analogous to the NLLs and 
NNLLs that are missing, in particular the one ~ 

Besides, it will be really worthwhile to use the recent calculations of the NNLO initial condition 
of the perturbative fragmentation function [27,28] and of the NNLO splitting functions in the non- 
singlet sector [29] to fully promote our formalism to NNLO/NNLL accuracy. This way we could 
explore whether the uncertainties on our predictions — especially the one due to the scale /%f ~ 
get finally reduced. If the theoretical uncertainties get substantially reduced, one may even think 
of extracting as(m 2 z ) from ^-fragmentation data, as already done from 5-decay spectra [9]. 

Our study could also be improved by including NNNLL contributions in the resummation 
exponents. Such corrections involve the coefficients B^ and D^: at present only B^ is 

exactly known. In particular, we expect a relevant effect due to the possible inclusion of D^ 3 \ the 
0(a%) coefficient of function D(as), which resums large-angle soft radiation in the initial condition. 
In fact, we have observed that the scale of as in the resummed initial condition is smaller than 
in the coefficient function, and therefore power corrections are more relevant. Also, our model 
enhances the coefficients of a s from the third order on. We can therefore employ our model 
to implement non-perturbative corrections to Drell-Yan or Deep Inelastic Scattering processes, 
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whose coefficient function was recently calculated to NNNLO accuracy [50]. In e + e~ annihilation, 
threshold contributions to the NNNLO coefficient function were computed in [51]; the full 
corrections have not been calculated yet. 

From the theoretical viewpoint, we plan to investigate in more detail the issue of the Mellin 
transform of the resummed cross section and the power corrections which are inherited by the B 
spectrum when one uses the effective as(k 2 ) and does the longitudinal-momentum integration in 
an exact or approximated way. Within our model, we chose to perform it exactly, driven by the 
results in [9,44], but nonetheless we believe that a thorough study of this point, along the lines 
of [35,36], should be necessary. 

Other issues which we plan to investigate in detail are the treatment of the higher orders of 
the effective as and the comparison between time- and space-like coupling constants. We have 
obtained reasonable agreement with the data by using the time-like effective coupling constant 
and taking the powers a's(k 2 ) after performing the integral of the dispersion relation. Nonetheless, 
this conclusion is somehow related to the level of approximation of the perturbative calculation 
which we have employed, and may not hold if we used a different accuracy. Therefore, a more solid 
understanding of these aspects of our model is mandatory. 

Furthermore, our model has some built-in universal features implying that, to put it into a 
stringent check, one should consider several observables from different processes. To this goal, 
however, one might need both NNLL resummation from the theoretical side and accurate data on 
the experimental one. Shape-variable data at the Z° peak, for example, are rather accurate, but 
for these quantities NNLL resummation is still in progress, hence preventing an analysis within our 
model. In [9] this model was applied to semi-inclusive B decays, where the situation is somehow 
complementary with respect to shape variables in e + e~ annihilation: NNLL threshold resummation 
is well established, but data are not very accurate yet, because of large backgrounds. 

Finally, we could also push our model to the 'low-energy' direction, to try to describe, for 
example, charm production at the Z° peak or below it, where accurate experimental data are 
available. Due to the large value of asijn 2 ) ~ 0.35 and to a soft scale S ~ m c (l — x), smaller by a 
factor of three than in 6-production, we expect a full NNLL/NNLO analysis to be necessary. The 
comparison with D-hadron data will be crucial to investigate possible deviations from our model. 
In fact, an extension of the formulation here presented may consist in adding to as{k 2 ) a correcting 
term: a' s (k 2 ) = as{k 2 ) + 5as{k 2 ). This way, as{k 2 ) will still be the effective coupling constant here 
discussed, while 5as(k 2 ) will depend on the hadronic state which we wish to describe, and vary, 
e.g, from mesons to baryons or from £?'s to L>'s. Investigating the charm sector may therefore help 
in shedding light on our model and understanding whether such a correcting term is necessary or 
not. This study is in progress as well. 
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A Numerical evaluation of Mellin transforms 



In this appendix we present a method for computing numerically the Mellin transform, through 
the Fast Fourier Transform (FFT), of a function of the form 



f(y) 



y 



(A.l) 



where y = 1 — z and ip(y) is a regular function of y or a function having at most a logarithmic 
singularity for y — > 0. 

To calculate the integer moments of the cross section <jn — typically the first few moments with 
A" = 2,3, 4, 5, • • • — the integral defining the Mellin transform does not present any convergence 
problem and can be done directly. In order to obtain the cross section in x-space, it is however 
necessary to perform an inverse Mellin transform by integrating along a vertical line in A^-space. 
The numerical computation of <jn in the complex A^-plane is non trivial for Im AT 1 because the 
kernel z N_1 develops fast oscillations with z, which affect the convergence of the integral. 



In detail, we have to compute the integral 



9n 



dy{l-y) N - l f{y) 



idy 
o y 



y) 



N-l 



<p{y)- 



for A" lying on a vertical line in the complex plane, 

N = c + iv, 



(A.2) 



(A.3) 



with c > and v real. It is convenient to treat analytically the infrared cancellation between real 
and virtual contributions related to the first and the second term in square brackets on the r.h.s. 
of Eq. ( 1A.2I) . We then take a derivative with respect to v: 



dg, 



c+iu 



dv 



Jo y 



(A.4) 



Note that the infrared singularity 1/y is now regulated by the ln(l — y) factor coming from the 
differentiation. 



In order to express 7 C (V) as the Fourier transform (FT) of some function, we change variables 



to 



t = - ln(l-y), 



(A.5) 
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and express Eq. (1A.4I) as an integral over t: 

lc {p) = / dte~ ivt ^ c {t), (A.6) 



with 

Mi) = ~i T^—t e ' Ct e i • ( A - 7 ) 

1 — e 1 1 J 

7 C is therefore the Fourier transform of the function ip c (t) defined above and vanishing for t < 0. 
For a fast numerical evaluation, it is convenient to transform the FT above into a Fast-Fourier- 
Transform (FFT). We cut the improper integral at a large but finite t mSuX : 

lc (v) ~ / dte~ lvt il) c {t). (A.8) 
Jo 

In practice, because of the exponential dependence y ~ exp[— t], we found that t max = 20 -j- 30 
already gives a good accuracy. The above integral can be easily approximated by means of a 
constant sampling: 

At = (A.9) 

n 

where n is the number of points in which the function tp c (t) is evaluated. We have found that an 
accuracy O(10 -3 -r 10~ 4 ) is reached already with n = 10 4 -j- 10 5 points. We then have: 



n— 1 r / -y 



k + -)At 



(A.10) 



7c (i/) ^ e~ t/2lyAi At ^e- ivfeA V c 

fc=0 

The Mellin transform is obtained by a numerical integration of 7 C (V): 

S'c+iv = 9c + [ dv'~i c {y'). (A.ll) 



g c is the Mellin transform in the fixed point c on the positive axis; it is a constant which is 
computed directly just once. Eq. (1A.11I) is our final result for the numerical computation of the 
Mellin transform. 

Let us now discuss the numerical computation of the inverse Mellin transform with a similar 
method: 

rc+ioc A AT 1 r+oo Ay 

M = / f L .{l-v)-' f 9N = / ^{l-yr v 9c +w . (A.12) 

Jc-ioo I'm (1 — y) c J-oo 2tt 

Since f(y) is real, one immediately obtains the following property of its Mellin transform: 

(9nT = 9n*. (A.13) 



The inverse Mellin transform can therefore also be written as: 

1 

7r(l -y) 
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1 r°° 

f(y) = -n T7 Re / due~ ivla ^g c+iv . (A.14) 



Apart from constant factors, Eq. ( 1A.14I) express the inverse Mellin transform as the inverse Fourier 
transform v — > — ln(l —y) of g c +iu- The transformation to the Inverse Fast Fourier Transform can 
be made as in the direct case. 

We have implemented the above algorithm within the Mathematica System with a gain of CPU 
time by over an order of magnitude with respect to the direct numerical evaluation. A typical 
run on a standard PC takes 0(1) minute. We have also checked our numerical results with direct 
integration in Fortran. 
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